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Advances in pure optical trapping techniques now allow the creation of degenerate Bose gases with 
internal degrees of freedom. Systems such as 8l Rb, 39 K or 23 Na in the F = 1 hyperfine state offer an 
ideal platform for studying the interplay of superfluidity and quantum magnetism. Motivated by the 
experimental developments, we study ground state phases of a two-component Bose gas loaded on 
an optical lattice. The system is described effectively by the Bose-Hubbard Hamiltonian with onsite 
and near neighbor spin-spin interactions. An important feature of our investigation is the inclusion 
of interconversion (spin flip) terms between the two species, which has been observed in optical 
lattice experiments. Using mean-field theory and quantum Monte Carlo simulations, we map out 
the phase diagram of the system. A rich variety of phases is identified, including antiferromagnetic 
(AF) Mott insulators, ferromagnetic and AF superfluids. 


I. INTRODUCTION 

The question of the interplay of superfluidity and 
internal bosonic degrees of freedom dates back many 
decades. At a purely conceptual level, the most 
straightforward issue is whether the two internal 
components move in, or out of phase. A generalization 
of Bogoliubov’s treatment to multiple species addressed 
this question and demonstrated that a neutral mixture 
of two species of charged bosons supports plasma- 
type excitations with oscillating charge density and also 
free-particle oscillations associated with mass density 
oscillations [1]. An extension of this work to finite 
temperatures considered dilute mixtures of (unstable) 
6 He in 4 He [2], Coupling between bosonic species 
was also shown to imply that superfluid motion of one 
component would result in a “drag effect” in which the 
second component is also set in motion [3]. 

Although different theoretical and experimental 
motivations were presented for studying multi- 
component bosons, in this early work, the physical 
system which was probably considered in the most detail 
was spin polarized hydrogen [4-6], where a large external 
magnetic field prevents recombination into molecules, 
and the smallness of the kinetic energy relative to the 
binding energy permits treatment with boson statistics. 
The key observation[7] was that these bosons reside 
in two low-lying hyperfine states, thus allowing for 
possible additional symmetry breaking associated with 
their relative occupation. The populations of the states 
were measured with electron spin resonance [8], and the 
nature of excitations of the “spin” degrees of freedom 
was shown to range from phonon-like, to free-particle- 
like with energy gaps, to resembling spin waves [9]. For 
a review, see Ref. [10]. 

Lattice models were also investigated. A mean field 


treatment[ll] of the hard-core limit of two component 
bosons focussed on the effect of “antiferromagnetic” 
interactions, i.e. a repulsion V between the different 
bosonic species on near-neighbor sites of a bipartite 
structure. Besides promoting an insulating phase where 
bosonic species alternate in a regular pattern, this 
interaction was found also to disrupt the “symmetrical 
condensate” [12] in which only a single superfluid species 
occurs, and allow for a superfluid phase in which both 
species are present. It was shown that, as a function of 
temperature T, two successive second order transitions 
can occur. For sufficiently large V, as T is lowered, 
the bosons first form a symmetric condensate and then, 
at a distinct, lower temperature, the asymmetrical 
condensate appears. We will show here that the soft-core 
lattice problem exhibits certain similarities with these 
hard core phase diagrams. 

Beginning in the late 1990’s, the properties of 
multicomponent boson systems became of renewed 
interest due to applications to ultracold quantum 
gases. Thanks to the all-optical trapping technique[13], 
hyperfine states of 87 Rb, 39 K or 23 Na in optical traps 
could now be used to realize interesting magnetic 
states[13-15]. In an optical lattice, it has been observed 
that atoms confined on the same lattice site exhibits 
collisions that could change their spin states[16]. Such 
a system can be effectively described by a multi- 
component Bose-Hubbard model with appropriate values 
of the intra- and inter-component interactions, and 
spin-conversion matrix elements[17, 18]. Due to the 
competition between spin species, the multi-component 
Bose-Hubbard model is expected to host novel phases[19- 
21] that are absent in the one-component Bose-Hubbard 
model[22]. 

Motivated by these experimental and theoretical 
developments for spinor bosons in optical lattices, we 
will study two component (spin-1/2) bosons on a 


two-dimensional lattice. We consider a very general 
Hamiltonian which includes not only on-site repulsion, 
but also near-neighbor interactions and interconversion 
between the species through a spin-spin coupling. We 
begin with a mean field theory (MFT) treatment which 
reveals a rich variety of magnetic patterns (unpolarized, 
ferromagnetic, and antiferromagnetic) accompanying the 
Mott and superfluid phases. Quantum Monte Carlo 
(QMC) calculations then are used to explore the phase 
diagram more exactly. In addition to showing that 
many aspects of the interplay between superfluidity and 
magnetism suggested by MFT persist, we also show that 
the order of the chemical potential driven superfluid- 
Mott phase transition depends on which Mott lobe is 
being considered, and even on whether commensurate 
density is being approached from above or below. 

High precision QMC work in two and three dimensions 
in the absence of interconversion has previously 
demonstrated the existence of different Mott and 
superfluid phases, distinguished by their patterns of 
charge and spin order [23]. The possibility of mixing 
‘heavy’ and ‘light’ bosonic species (‘mass imbalance’) 
introduced additional phenomena like ferromagnetic, 
phase separated states, and ‘entropy squeezing’, in which 
the heavy species is in a Mott phase while the light 
species is superfluid and can act as a heat reservoir to 
absorb entropy [24]. The effects of interconversion on 
these phenomena is one of the topics of the present work. 

While we will explore here the phase diagrams for quite 
general values of the kinetic and interaction energies, 
we note that the precise quantitative form of the 
effective (pseudo) spin interaction potential for ultracold 
bosonic and fermionic atoms can be computed using 
the “degenerate internal state approximation” [25] . At 
a basic conceptual level, the coupling is similar to a true 
spin interaction in which the magnetic field produced 
by one spin couples to the second spin, but there are 
important differences. One of these is that, because 
the hyperfine states are not “real spin”, they are not 
generators of rotations, and hence there is no reason to 
expect an isotropic (“Heisenberg-like”) form J Si • £> 2 . 
Instead, the energy can be Ising or XY in character, 
and indeed the precise form depends on the scattering 
lengths of binary atom-atom collisions in the presence of 
an external field. 


II. THE SPIN-1/2 MODEL WITH 
NEAR-NEIGHBOR SPIN INTERACTIONS 


Here we are interested in the spin-1/2 Bose-Hubbard 
Hamiltonian[17, 18] with the spin interactions extended 
to near-neighbor sites 


H = ~tj2 (EE + EE) -^E f Ai + ^ E ^ 
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In the above equation, b' ai (b ai ) creates (annihilates) a 
pseudo-spin a boson on site i of an L x L square 

lattice under periodic boundary conditions, h\ = n-j-i + 
hi i, and 5“ ( a = x,y,z) is the spin operator defined as 

srEE^C'i i. (2) 
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where < 7 “^, are the Pauli matrices. The parameters t 
and y correspond to the near-neighbor (NN) hopping 
amplitude and chemical potential respectively. We use 
t = 1 as the unit of energy. Uq is the contact interaction, 
while U 2 and V are on-site and NN spin-spin interactions 
respectively. 

Using the representation Eq. (2) and taking the NN 
spin-spin interaction to be along the 2 -axis only [26], we 
arrive at the following model Hamiltonian 
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It can be seen that the onsite spin coupling U 2 has 
two roles. First, it shifts the strength of the contact 
interaction between opposite spins Second, U 2 is 

also the matrix element of the conversion process which 
turns two identical bosons into the opposite spin species 
when they meet at the same site. 

Eq. (3) is the Hamiltonian that will be studied in this 
work. It will be solved using MFT and exact stochastic 
Green function (SGF) quantum Monte Carlo technique. 
While the SGF QMC method can treat the contact spin 
interactions or the NN spin-spin couplings separately, a 
sign problem would arise if both terms were retained 
due to the presence of interconversion matrix elements in 
both. For this technical reason, we drop the conversion 
matrix elements of V terms in Eq. (1) and study Eq. (3). 
The retention of the 2 -axis term gives important insights 
into the effects of the NN spin-spin interactions, in the 
same spirit that the t-J z Hamiltonian provides initial 
clues into the more general rotationally invariant t-J 
model. 


3 


We focus on the case where U 2 > 0 and V < 0, 
i.e., antiferromagnetic NN spin couplings. In general, 
the value of U 0 , U 2 and V will depend on details of 
the system (for example, scattering length between the 
atoms, polarization of the laser waves forming the optical 
lattice, detuning from the internal atomic transition 
etc. [17]). Here we adopt the parameter regime studied 
in Ref. 17 where, based on known values of s, Rb 
and 23 Na scattering lengths and on laser wavelengths 
corresponding to the D\ resonance, U 2 is typically an 
order of magnitude or more smaller than Ho- In the 
current work, we take the value U 2 /U 0 = 0.1. For the 
NN coupling V, because interaction strength typically 
decreases with distance, we assume the value |U|/f7o < 
0 . 1 . 

When V = 0 in Eq. (3), the system was studied 
extensively by MFT[17, 18] and QMC methods in one 
and two dimensions [27, 28] . In 2D and U 2 > 0[29] , the 
ground state of the Hamiltonian features three phases: a 
ferromagnetic superfluid (FMSF), an unpolarized Mott 
insulator (MI) at even commensurate densities, and a 
ferromagnetic Mott phase at odd commensurate Filings. 
For negative U 2 , it was found that the ground state never 
polarizes[27, 28]. 


III. MEAN FIELD THEORY 

A. Decoupling Mean Field theory 

The mean-field scheme employed in the present work 
is developed in Ref. 30 and 31. The method is based on 
rewriting the Hamiltonian as a sum over local terms that 
can be solved exactly for a fixed number of bosons. To 
incorporate the hopping terms, one introduces uniform 
SF order parameters (b'^) = (S CTi ) = ip G . Since we are 
interested in equilibrium states, the order parameters ip a 
can be chosen to be real. Using this ansatz, the kinetic 
energy terms, which are non-diagonal in boson creation 
and destruction operators, are decoupled as 

~(kti + &aj)VV + V£, (4) 

where in the last line we have dropped the terms that 
have products of fluctuations in bosonic operators on 
different sites. 

To treat the NN spin interactions in the same 
decoupling scheme, we decompose the square lattice into 
two disjoint sublattices A and B and introduce real 
magnetic order parameters (S A ), (S B ) on sublattice A 
and B respectively. Under the MF approximation, the 


spin-spin interaction term now becomes 

\V\J2SpS? * 2 c |Vj ( s?(s z B ) - U& A )(fr B )) 

(ij) V J 

(5) 

As before, in this form, terms that describe products 
of fluctuations in S Z A and S B are ignored. z c = 4 is 
the coordination number of the square lattice. We have 
also assumed that the magnetic order parameter on each 
sublattice is uniform. With these approximations, the 
final MF Hamiltonian becomes two coupled local ones 
for sublattices A and B: 

H e = - Zcty^XKe. + i>ae)ipae + z c t ^ ipatipa-e 

a a 

+ ~^h n {h n - 1 ) + ^hieinie - 1 ) 

+ ~2~ + ^- c -) + (^0 + U 2 )nnhn 

- Knu + fi U ) + z c \V\St(SI) - ^M(5|)(Sf), 

( 6 ) 

where l = A, B (l = B, A), and a =t>4- (d = 
|,t)- The coupled Hamiltonians are solved at zero 
temperature by iteration. Starting with an initial guess 
of order parameters ip a t and (5|), Ha and H B can be 
diagonalized numerically within the bosonic occupation 
number basis {|n^,nj,f}} truncated at n a i < fV max . 
Order parameters are then updated with respect to 
the new MF ground state. This procedure is repeated 
until ip„t, (Sp and the ground state energy E g are 
converged. We typically choose lV max = 14 to ensure 
that convergence is independent of -/V max . Multiple initial 
configurations are also used to verify that the converged 
MF solution do not depend on initial conditions. We 
benchmark our MF program by computing the phase 
diagram of Eq. (6) with V = 0. The results are in 
agreement with previously published data[27]. 

Different MF phases are classified by the corresponding 
order parameters. For example, a superfluid is 
characterized by finite total superfluid density 

Ps ,1 — PsM + Ps,U = + ipfi- ( 7 ) 

The Mott insulator, on the other hand, is defined by 
zero superfluid density p s / = 0 and zero compressibility 
dpc/dp = 0, where 

Pi = P^t + Pie- (8) 

To examine magnetic order, we compute the expectation 
value of S z with respect to the converged MF solution. 
In a Mott phase, this is 

S z e = \{n\e~fiu), (9) 
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FIG. 1. (Color online) Schematic diagram of ferromagnetic, 
antiferromagnetic, and unpolarized states at different 
commensurate fillings. Blue (+) and red (—) circles represent 
spin-up and spin-down components respectively. 


where h a n is the density operator. In principle one can 
use Eq. (9) in the SF phase, and the conclusion should 
remain the same. Here we follow the convention in 
Ref. 17 and 32 and compute the magnetization in the 
SF phase defined as 


qz _ 


1 t ~ 

2 p s ,e 


( 10 ) 


which merely measures the SF population difference 
between the two spin components. 

Figure 1 summarizes schematically the possible 
magnetic structures at three commensurate fillings. For 
example, the state (SF or MI) is ferromagnetic (FM) if 
one of the spin components dominates the population 
throughout the lattice. An unpolarized state has 
both spin components equally occupied on every lattice 
site. An antiferromagnetic (AF) state is realized when 
sublattices A and B are dominated by different spin 
species. 


■o Pa ' Pb -O Ps,a * Ps,b 




FIG. 2. (Color online) Features of the MF ground state as a 
function of chemical potential p/Uo at t/Uo = 0.02, U 2 /U 0 = 
0.1 and \V\/Uo = 0.02. Total particle and total superfluid 
densities on l = A and B sublattices are plotted in the top 
panel. Magnetic order parameters Sf and S/ t are shown in 
the bottom panel. In both figures, the vertical dashed lines 
divide the phases into various zones based on their magnetic 
structures labeled by the roman numerals: I - unpolarized, 
II - AF, and III - FM. Note that the changes in sign of S\ 
(Sg) in the first Mott plateau are a trivial reflection of the 
degeneracy of the two MF solutions. 


B. Mean Field Results 

Properties of the MF ground state are shown in Fig. 2 
for t/Uo = 0.02, U 2 /U 0 = 0.1, and \V\/U 0 = 0.02. Total 
particle and SF densities are plotted in the upper figure 
as functions of p/Uo- The density develops three well 
defined plateaux at p = 1, 2, and 3. These plateaux 
correspond to Mis because the compressibility dp/dp = 0 
and the SF density also vanishes. The SF phase resides in 
between the Mott insulators. It can be seen that pt and 
p s j change discontinuously when one enters and leaves 
the first Mott plateau. This is the signature of a first 


order phase transition. Likewise, the transition is also 
first order as one enters the p = 3 MI from below. On 
the other hand, pc and p s ^ change continuously on both 
sides of the second Mott plateau, indicating that the MI- 
SF transition is second order. 

The lower panel of Fig. 2 summarizes MF magnetic 
structures for t/Uo = 0.02, U 2 /U 0 = 0.1, and \V\/Uq = 
0.02. Within the p = 1,3, and 4 plateaux, the magnetic 
order parameter on sublattices A and B are equal but 
have opposite signs S\ = — S%. This shows that these 
Mis are antiferromagnetic. In contrast, the second Mott 
insulating region is non-magnetic. In the SF region, 
magnetic properties are plotted by blue symbols (dots 
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and empty circles). The SF between the p = 1 and 
p = 2 Mott regions has two different magnetic natures: 
unpolarized and fully polarized. The transition between 
them is first order. This is also indicated in the upper 
panel by a discontinuity in p s ^. Most interestingly, the 
SF above the p = 3 Mott region shows antiferromagnetic 
structure. 

By carrying out the self-consistent MF calculation at 
different t/Uo (or p/Uo) values, the p-t phase diagram 
can be constructed. Results for \V\/Uq = 0.02 and 
0.08 are plotted in Fig. 3. Here red (solid) and blue 
(dashed) curves represent first and second order phase 
transitions respectively. Comparing with the V = 0 
phase diagram[27], there are several notable changes due 
to the presence of NN spin-spin couplings. 

The magnetic structure of the first and third 
Mott lobes changes from being ferromagnetic to 
antiferromagnetic. At p = 1 or 3, one of the spin 
components dominates the population. As a result, the 
MF ground state energy can be lowered by forming an AF 
pattern. At p = 2, on the other hand, the onsite coupling 
term can be effectively avoided by equally populating 
both spin species on every site if \V\/Uo = 0.02 is small. 
By raising \V\/Uo to 0.08, the second lobe also becomes 
antiferromagnetic. This is because the energy gained 
by forming an AF state compensates the energy cost of 
onsite coupling terms at large \V\/Uo values. 

When V = 0, the MI-SF phase transition is continuous 
except for the tip of the second Mott lobe. The transition 
is known to be first order for 0 < U^/Uq < 0.25. [27] 
Here we find that the transition becomes first order 
due to the change of magnetic property in the p = 1 
and 3 Mott lobes. Note that above the third lobe, the 
antiferromagnetic Mott insulator to antiferromagnetic 
superfluid (AFSF) transition remains continuous. At 
\V\/U 0 = 0.08, the bottom half of the phase boundary 
enclosing the Mott insulators is first order; while the 
upper half becomes continuous. 

Regarding the magnetic structure of the SF phase, 
it was found that the SF is always polarized if V = 
0[27] . With the presence of NN spin-spin couplings, an 
unpolarized SF emerges near the Mott lobes, particularly 
at small t/Uo values. An exception to this observation 
is found at p > 3 where an AFSF phase occupies the 
region between the third and fourth Mis. At \V\/Uq = 
0.08, the AFSF region expands dramatically to large 
hopping regions, and to chemical potential values as low 
as p/Uq ~ 0.6. This AFSF is a supersolid phase since 
it exhibits simultaneous diagonal and off-diagonal long 
range order. 

Recall that in the original Bose-Hubbard model [22] or 
in the case V = 0 in Eq. (6) [27, 28], the SF phase extends 
all the way to t/Uo = 0. Fig. 3 shows that this is no longer 
the case when V is turned on. The system undergoes a 
series of first-order transition between Mis at small t/Uo- 



FIG. 3. (Color online) MF phase diagrams obtained by 
solving the coupled Hamiltonian Eq. (6) for U 2 /U 0 = 0.1. The 
NN spin interaction strength \V\/Uo is (a) 0.02 and (b) 0.08. 
First order and continuous phase transitions are represented 
by red (solid) and blue (dashed) curves respectively. “AF”, 
“FM”, and “unpolar” indicate magnetic structures. At 
\V\/U 0 = 0.02, an AFSF phase emerges between the p = 3 
and 4 Mott lobes. With increasing \V\/Uo, the AFSF region 
expands significantly. Moreover, the p = 2 lobe changes 
its nature from being non-magnetic at \V\/Uo = 0.02 to 
antiferromagnetic at \V\/Uo = 0.08. In panel (a), the vertical 
dashed line indicates the location of t/Uo where Fig. 2 is 
plotted. 


IV. EXACT QUANTUM MONTE CARLO 
STUDY 

In this section, we solve the model Eq. (3) exactly on 
finite lattices by using Stochastic Green Function (SGF) 
QMC[33]. The SGF method is a finite-temperature 
continuous time QMC technique that can be formulated 
in either the canonical or grand canonical ensembles. 
The SGF algorithm can solve a large class of lattice 
Hamiltonians that can be written as H = V — T, where V 
is diagonal in the Fock basis (subject to the model type) 
and T has only positive elements[33]. The technique has 
also been applied to the V = 0 case of Eq. (3) in one and 
two dimensions [27, 28]. 

In our simulations, the temperature is set at f3t = 
2 L for a lattice with linear dimension L. The chosen 
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FIG. 4. (Color online) Particle density p (square and triangle) 
and SF density p s (diamond) versus chemical potential 
computed using the SGF QMC method for t/Uo = 0.02. 
For comparison, we have implemented canonical and grand 
canonical ensembles in the simulations. Both data sets 
are acquired on a L = 6 lattice at temperature ft = 12. 
Conversion term and NN spin interactions are U 2 /U 0 = 0.1 
and \V\/U 0 = 0.02 respectively. The canonical ensemble 
results are shifted upward by 0.5 for clarity. The signature of 
first order transition (discontinuous jump in p and p s in the 
grand canonical ensemble curve, and negative compressibility 
dp/dp < 0 in the canonical ensemble data) can be seen at 
both ends of the first Mott lobe and the bottom boundary of 
the third lobe. 

temperature is typically low enough to ensure that the 
results are converged to the ground state limit. In 
some cases, we select j3t = 4 L to reach convergence. 
We benchmarked the SGF algorithm by comparing with 
exact diagonalization data for a small cluster. The SGF 
and exact results are in agreement within statistical 
errors. 


A. Phase diagram 



FIG. 5. (Color online) Exact phase diagram of the 
Hamiltonian Eq. (3) obtained using the SGF QMC technique 
on L x L lattices at temperature f) = 2L. The onsite 
spin coupling is fixed at U 2 /U 0 = 0.1. The NN spin- 
spin interaction \V\/Uo is set at (a) 0.02 and (b) 0.08. 
Statistical uncertainties are typically smaller than the symbol 
size. The vertical dashed line is the location where Fig. 4 
is plotted. Mean-field phase boundaries are also presented 
(dashed curves) for comparison. In (a), the vertical heavy 
line in the p = 1 Mott lobe near the tip is the location of 
magnetic transition (c.f. Fig. 7). In (b), AFSF is predicted 
to exist in the region above the p = 3 MI. 


To construct the exact phase diagram, we compute the 
total particle density p and SF density p s as functions 
of chemical potential. In canonical ensemble SGF 
simulations, the total particle number N is fixed, and 
we derive the chemical potential via 

p(N) =E(N + 1)-E(N), (11) 

where E(N) is the total energy of N bosons on an L x L 
lattice. To access the SF density, we use the formula 
proposed by Pollock and Ceperley[34], which relates p s to 
the winding number W. However, due to the conversion 
term in the Hamiltonian Eq. (3), the numbers of spin 
t and l bosons are not conserved individually. As a 
consequence, the relevant winding number should take 
into account both spin components [35] and p s is given 


by the following formula 


+ wtf) 

2 dt/3L d ~ 2 


( 12 ) 


where d = 2 is the dimensionality, t is the hopping 
amplitude, and W-j- and Wj, are the winding numbers of 
spin f and j, bosons respectively. 

Figure 4 shows QMC results for p and p s versus p/Uo 
on the L = 6 lattice with t/Uo = 0.02, U 2 /U 0 = 
0.1, and \V\/U 0 = 0.02. We compare total densities 
p(p) measured using both grand canonical (square) 
and canonical (triangle) ensembles. Three plateaux 
can be observed at commensurate fillings. Since the 
compressibility dp/dp and superfluid density p s vanish 
in the plateaux, these regions represent Mott insulators. 
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In between the Mott insulators there is a SF with p s ^ 0. 
The agreement of the data for different ensembles acts 
both as a check of our codes and also as an assessment of 
finite size effects, since equivalence is expected only for 
sufficiently large lattices. 

In Fig. 4, the grand canonical ensemble particle density 
has a discontinuous jump when one enters or leaves 
the first Mott region. These jumps show that the MI- 
SF transition is first order. This is confirmed by the 
canonical ensemble data which clearly indicates negative 
compressibility dp/dp, < 0 in the same region. At the 
same time, the SF density shows a discontinuous jump. 
Likewise, the MI-SF transition near p/Uq ~ 2.1 is first 
order. The transition at p = 2 is second order as 
both quantities p and p s change continuously (within the 
resolution of our p/Uq grid) as a function of the chemical 
potential. MF predictions at t/Uo = 0.02 (cf. Fig. 2) are 
consistent with these exact QMC results. 

The QMC phase diagram is shown in Fig. 5 for (a) 
\V\/Uq = 0.02 and (b) \V\/Uq = 0.08. The onsite spin 
coupling strength is U 2 /U 0 = 0.1 in both cases. The 
corresponding MF phase boundaries (dashed curves) are 
also plotted for comparison. The QMC data are shown 
for p/Uq < 2.75 as it becomes increasingly difficult to 
reduce statistical errors for simulation at large chemical 
potential values. System sizes L = 6, 8, 10, and 12 are 
used, with little variation evident on the QMC phase 
boundaries. Overall, the QMC and MF phase boundaries 
are in good agreement, especially at small t/U q where the 
MF assumption works well. The deviation between the 
two approaches increases as one moves toward the tips 
of the Mott lobes where quantum fluctuations are large. 
Interestingly, at /Uq = 0.02 and near t/U 0 ~ 0.045, our 
QMC data reveal a magnetic phase transition inside the 
first MI. The transition is indicated by a thick black line 
in Fig. 5 and will be discussed in the next subsection. 
At \V\/Uq = 0.08, the QMC Mott insulating regions 
expand, which is consistent qualitatively with MF results. 
At small t/U 0 values, our QMC data also suggest the 
existence of a direct first-order MI-MI transition at both 
\V\/Uq = 0.02 and 0.08, confirming the MF predictions. 


B. Magnetic properties of the Mott lobes 

To study magnetic properties of the model in 
QMC simulations, we measure the real-space spin-spin 
correlation function along the 2 -axis 

Cs{r)=^J2^r+r'S; l ). (13) 

r' 

Figure 6 shows the results obtained for the L = 12 lattice 
at p = 1 (upper panel) and p = 2 (lower panel) with 
U 2 /Uq = 0.1, \V\/U a = 0.02, and t/U 0 = 0.02, i.e. 
inside the MI phases. The staggered correlation pattern 
displayed in Fig. 6 (a) shows that the first Mott lobe is 
antiferromagnetic. Similar results are also obtained for 


(a) p = l.O 
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FIG. 6. (Color online) Spin-spin correlation function C zz (r) 
for the (a) p = 1 and (b) p = 2 Mott insulators. The data are 
measured on an L = 12 lattice for t/Uo = 0.02, U 2 /U 0 = 0.1, 
and \V\/Uq = 0.02. The data show that the first lobe has AF 
order while the second lobe is non-magnetic. 


the third lobe. As indicated by Fig. 6 (b), the second 
Mott lobe at \V\/Uq = 0.02 is non-magnetic since only 
short-ranged correlation exists. 

In order to confirm that the p = 1 and 3 Mis have long- 
range magnetic order, we have also studied the scaling of 
the spin structure factor at (7r, 7 r) 

S*°(n,n) =£(-l) p C7 8 (r). (14) 

r 

If the state has a long-range AF order, then S zz (tt,tt) 
should scale as L 2 [36]. The results for the first Mott lobe 
at U 2 /Uq = 0.1 and \V\/Uo = 0.02 are plotted in Fig. 7 
(a) as a function of t/Uo. In this figure, the filled symbols 
represent S zz (Tr,n)/L 2 computed for L = 6, 8, 10, and 
12. The data confirm that the first Mott lobe has long- 
range AF order. By carrying out similar scaling studies, 
we have verified that the third lobe at U 2 /Uq = 0.1, 
\V\/Uq = 0.02 (cf. Fig. 7 (b)), and the p = 1, 2, 
and 3 MI phases at U 2 /Uq = 0.1, \V\/Uq = 0.08 also 
have long-range AF order. These findings confirm the 
MF predictions regarding the magnetic structure of Mott 
insulators at commensurate fillings in the parameter 
ranges studied. 
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FIG. 7. (Color online) Normalized antiferromagnetic spin 
structure factor S zz (n, 7 t)/L 2 and SF density p a as functions 
of hopping amplitude t/Uo- Here U 2 /U 0 = 0.1 and \V\/Uo = 
0.02. I 11 panel (a), the data are plotted near the tip of the 
p = 1 Mott lobe of the phase diagram Fig. 5. AF order 
is destroyed at t/Uo ~ 0.045, before the system becomes a 
SF near t/Uo ~ 0.06. The shaded region denotes the non¬ 
magnetic Mott insulator. In panel (b), the same observables 
are plotted at p = 3. Here the insulator-to-SF and magnetic 
transitions take place at the same critical values t/Uo ~ 0.39 
and is first order. 


p s = p g ^ 4 - The system is a L = 6 lattice with U 2 /U 0 = 
0.1, \V\/Uo = 0 . 02 , and /? = 12 . In case (a) where t/Uo = 0.1, 
the probability P(p SjCT ) centers at p s / 2, indicating that both 
spin components are equally populated in the SF state. In 
case (b) t/Uo = 0.5, there is an asymmetry between P(ps,t) 
and P(p s , 4 .) at p s = 3. This implies that the SF state has a 
finite polarization. 


C. Magnetic properties of the SF phase 


In addition to the AF order parameter S zz (tt,tt), we 
also show in Fig. 7 the total SF density as a function of 
t/Uo for \V\/Uo = 0.02. It can be seen that at p = 1, 
the SF density p s rises and becomes size-independent 
(indicating a true SF phase) at t/Uo ~ 0.06, a value that 
is consistent with the one found in Fig. 5 (a). Therefore, 
as one scans through t/Uo, the P = 1 Mott insulator 
undergoes a first order (indicated by the discontinuous 
jump in S zz (tt,tt)/L 2 ) magnetic phase transition at 
t/Uo ~ 0.045 before it becomes a SF. This magnetic 
phase transition is not captured by the MF theory. 

Figure 7(b) shows a similar analysis near the tip of the 
third MI phase at \V\/Uq = 0.02. It is found that the 
MI-SF transition is first order and takes place at t/Uo ~ 
0.037. However, no intermediate phase exists. 


Next we turn our attention to magnetic properties 
of the SF phase. MFT predicts three different types 
of SF: a FMSF, an AFSF, and an unpolarized SF. At 
\V\/U 0 = 0.02, the FMSF dominates the phase diagram. 
At a stronger NN spin coupling \V\/U 0 = 0.08, the AFSF 
becomes the major component(cf. Fig. 3). 

To verify these MF predictions, we first compute the 
SF density histogram P(p Sj(T ) for spin a =f and j. 
bosons. As shown in previous results[27, 28], P(p S)(T ) for 
a =f and j, are identical if both spin species are equally 
populated. On the other hand, P{p s> - (•) and P{p 8 , 4 .) would 
peak at different values of p Sy(7 if the superfluid develops 
polarization. 

At U 2 /U 0 = 0.1 and |Vj/[/o = 0.02, examples of the 
histogram are plotted in Fig. 8(a) for the L = 6 lattice at 
three commensurate densities p s = p s ,t + Ps,i = 1, 2, 3 
at t/Uo = 0.1, i.e. deep inside the SF phase. The figure 
shows that P(p SjCT ) are identical for both spin components 
at a given density and peaks at p s /2 , indicating no spin 
polarization. We find no FMSF in the parameter range 
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FIG. 9. (Color online) Upper panel : Antiferromagnetic 
spin structure factor S zz (n, n) and SF density p s versus t/Uo 
measured at p = 3.5. The onsite and NN spin interactions 
are U 2 /U 0 = 0.1 and \V\/Uq = 0.08 respectively. The system 
is a SF for t/Uo > 0.02. Within the SF phase, S zz (n,n) 
scales as L 2 . These results suggest that the SF has long- 
range AF order. Lower panel: Spin-spin correlation function 
C's(r) computed on the L = 12 lattice at the t/Uo = 0.0588, 
indicated by the dashed vertical line in the upper panel. 


shown in the \V\/Uq = 0.02 QMC phase diagram (top 
panel of Fig. 5). 

In order to search for the FMSF further, we carry 
out the simulation at much higher t/Uo values. One 
representative result of P(p St<T ) at t/Uo = 0.5 is depicted 
in Fig. 8(b) for the L = 6 lattice with U 2 /U 0 = 0.1 and 
\V\/U 0 = 0.02. The figure shows that at p s = 1, 2, and 
3, the histogram peaks at different locations for different 
spin species. These results suggest the existence of a 
spin polarized SF phase, albeit at a much higher hopping 
range than the MF prediction. A similar conclusion is 
reached at \V\/U 0 = 0.08 for the FMSF phase. 

To probe the AFSF phase, we compute the spin 
correlation function Eq. (13) in the superfluid phase. 
At |E|/17o = 0.02, results only indicate short-range AF 
correlations, and the corresponding scaling study of AF 
spin structure factor does not support any long-range AF 





FIG. 10. (Color online) Normalized spin structure factor 
at ( 71 ", 7r) (top panel) and superfluid density (bottom panel) 
versus chemical potential p/Uo computed on the L = 6 lattice 
with U 2 /U 0 = 0.1, \V\/Uo = 0.08. In the top panel, there is 
a sudden increase in S zz (tt,tv) near p/Uo ~ 2.2, 3.1, and 3.9 
for t/U 0 = 0.071. At roughly the same p/Uo values and the 
same t/Uo, a reduction in p s can been observed. 


order in the superfluid. 

We carry out the same analysis for |U|/t/o = 0.08, 
and the results at p = 3.5 are summarized in Fig. 9. 
The upper panel of the figure shows S zz (tt, 7 t)/L 2 as well 
as total SF density p s in a range of t/Uo values. The 
superfluid density data indicate the onset of superfluidity 
is at t/Uo ~ 0.02. In the SF phase, the AF spin structure 
factor is finite and scales as L 2 before it drops to zero at 
t/Uo ~ 0.08. These data combined therefore confirm the 
existence of an AFSF phase at \ V\/Uq = 0.08. This phase 
can be considered a supersolid phase since it exhibits 
similtaneous digaonal and off-diagonal long range order. 
In the lower panel of Fig. 9, we show a real-space spin 
correlation function result acquired on an L = 12 lattice 
at p = 3.5, U 2 /U 0 = 0.1, \V\/U 0 = 0.08, and t/U 0 = 
0.0588 (indicated by the vertical dashed line in the upper 
panel of Fig. 9). The staggered correlation function 
pattern demonstrates the long range antiferromagnetic 
structure of the SF phase. This long range AF order in 
the SF phase appeared as the NN repulsion was increased 
from \V\/Uq = 0.02 to \V\/Uq = 0.08. We have not, 
however, determined the value of \V\/Uo at which the 
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AFSF first appears. 

Figure 9 also indicates that at some values of t/Uo, 
the AF order vanishes and the SF becomes a normal 
superfluid. To estimate the exact phase boundary of 
this AFSF to normal SF transition, we have conducted a 
series of grand canonical and canonical SGF simulations 
and extracted S zz (tt 1 tv) as a function of p/U 0 for several 
t/Uo values. A set of data is presented in Fig. 10 for 
the L = 6 lattice with U 2 /U 0 = 0.1, \V\/Uq = 0.08, 
and 2 < p/Uo < 4. It can be seen from the figure that 
the transition from a normal SF to AFSF takes place at a 
chemical potential value much higher than the MF result, 
and the critical p/Uo increases with t/Uo. We have done 
similar calculations for other system sizes. The estimated 
AFSF phase boundary is shown in Fig. 5 (b). 

Finally, we would like to remark that as one scans the 
chemical potential in the range 2 < h/Uq < 4 at t/Uo = 
0.071, a reduction in the SF density at p/Uo ~ 2.2, 3.1 
and 3.9 can be observed in Fig. 10. Correspondingly, 
S zz (rv, tv) rises rapidly near these p/Uo values. Because 
t/Uo = 0.071 is just outside the tip of the third Mott 
lobe, the reduction in p s at p/Uo = 2.2 is caused by the 
proximity effect of the third MI phase. The reduction in 
p s at p/Uo ~ 3.1 indicates indirectly the location of the 
tip of the fourth AF Mott lobe (and potentially the fifth 
at p/Uo ~ 3.9). 

V. CONCLUSION 

In this work, we have used the site-decoupling MF 
theory and the exact SGF QMC algorithm to study the 
ground state phase diagram of the Bose-Hubbard model 
with onsite and NN spin-spin couplings Eq. (3). The 
SGF approach allows us to treat terms which interconvert 
the two bosonic species. Previous study[28] have shown 
that the Hamiltonian at V = 0 and positive U 2 has 
three phases: a ferromagnetic Mott insulator at p = 1 
(and all odd Mott lobes), an unpolarized Mott phase 
for p = 2 (and all even commensurate densities), and 


a ferromagnetic superfluid. 

In the presence of NN interactions \V\ ]C(ij) 5'?, the 
magnetic structures found at V = 0 are profoundly 
modified. In particular, at \V\/Uo = 0.02, the Mott 
lobes at p = 1, 3 become antiferromagnetic. The Mott 
phase at p = 2 is a spin-singlet state. The superfluid 
phase becomes unpolarized for t/Uo 1$ 0.5. By increasing 
the strength of NN spin coupling to \V\/Uo = 0.08, the 
second Mott lobe also becomes antiferromagnetic, and, 
most interestingly, an AFSF (a supersolid phase) emerges 
at high fillings. 

At the \V\/Uo values studied, the MF and exact QMC 
results are in good agreement, particularly at small t/Uo 
values (deep inside the Mott phase) where quantum 
fluctuations are small. Moreover, the site-decoupling 
MFT is able to capture correctly the magnetic structure 
of the Mott insulators and predict the existence of AFSF. 
The order of MI-SF phase transition is also verified by 
the exact results. 

Just as initial qualitative studies of the single species 
boson-Hubbard model were followed by quantitative 
comparisons with experiment [37, 38], a natural next 
step here will be to do similar modeling of multi- 
component bosonic optical lattice experiments. However, 
the complication introduced by the effect of a trap, 
which in the single species case manifests itself as the 
coexistence of superfluid, Mott insulator, and normal 
phases as p, U/t and T/t vary across the cloud, will be 
even more challenging, since the possibility of magnetic 
order introduces additional phases which might coexist 
in the presence of a confining potential. 
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